Changes in Alpine Butterfly Communities during the Last 40 Years

Simple Summary The Alps are among the most vulnerable ecosystems to climate change, as these modifications take place at a faster pace at higher elevations. Since butterflies are ideal model systems to investigate species responses to climate and habitat changes, we monitored a butterfly community of the Valasco valley in the Maritime Alps (NW Italy) by three sampling events covering a period longer than 40 years. We observed an overall increase in mobile, tolerant and thermophilous species, which eventually might cause an overall loss of community distinctiveness. The variations observed in the butterfly community can be explained by the notable increase in maximum temperatures and the reduction of grasslands, along with the increase of woodlands, supporting the hypothesis that local warming and land-use changes have ultimately affected the butterfly community composition. Abstract Our work aims to assess how butterfly communities in the Italian Maritime Alps changed over the past 40 years, in parallel with altitudinal shifts occurring in plant communities. In 2019, we sampled butterflies at 7 grassland sites, between 1300–1900 m, previously investigated in 2009 and 1978, by semi-quantitative linear transects. Fine-scale temperature and precipitation data elaborated by optimal interpolation techniques were used to quantify climate changes. The changes in the vegetation cover and main habitat alterations were assessed by inspection of aerial photographs (1978–2018/1978–2006–2015). The vegetation structure showed a marked decrease of grassland habitats and an increase of woods (1978–2009). Plant physiognomy has remained stable in recent years (2009–2019) with some local exceptions due to geomorphic disturbance. We observed butterfly ‘species substitution’ indicating a general loss in the more specialised and a general gain in more tolerant elements. We did not observe any decrease in species richness, but rather a change in guild compositions, with (i) an overall increased abundance in some widespread and common lowland species and (ii) the disappearance (or strong decrease) of some alpine (high elevation) species, so that ‘resilience’ could be just delusive. Changes in butterfly community composition were consistent with predicted impacts of local warming.


Introduction
Climate studies strongly suggest that atmospheric alteration is already occurring and that changes in atmospheric composition are altering weather and climate processes. more promptly than many other organisms to environmental change [47]. Since climate change can alter species dynamics, results obtained by the study of changing butterfly communities in a varying climatic scenario will provide a starting point for research on a wide variety of other herbivorous insects [48][49][50]. Butterflies, as well as biotic components in general, may adapt to changing environments showing phenological plasticity and/or genetic changes (see [51] for a review on invertebrates), as well as by shifting their spatial distribution moving towards favourable conditions [52]. Thus, as the climate warms, adults tend to emerge earlier [53,54], flight periods last longer [54] and numbers of generations in multivoltine species increase [55]. In mountain ecosystems, upwards shifts of species linked to colder conditions can be observed [8,52]. Thermophilous species could profit from global warming and extend their occurrence by spreading northwards or to higher elevations. However, species movements occur only if suitable and relatively uninterrupted habitats are available (mainly represented by open areas in the Alps), and often mirror range contractions [56][57][58]. In the worst cases, global warming and loss of suitable habitat drive local extinctions, but species responses might also lead to less severe changes at the community level [59]. In a long-time series study (1840-2013), Habel and colleagues [59] found an overall impoverishment of grassland butterfly communities in southern Germany, with a decrease in the proportion of specialists (already endangered species) and lowdispersal butterflies. Even considering a short time frame, Cerrato et al. [60] observed quite significant variations in butterfly assemblages of the NW Italian Alps which can be partially explained by both global warming and microhabitat changes. These authors found an overall increase of generalist and highly vagile species in five years, leading to an incipient community homogenisation [60].
Although, in the mountains, climate alterations can lead to rapid changes in ecosystems and in the way species occur in space and interact with each other in communities, long-term studies on variations of invertebrate assemblages are surprisingly rare (but see [60]) in the Alps. Therefore, in 1978, 2009 and 2019 we sampled butterfly species occurring in seven plots of the SW Italian Alps, located at elevations ranging from 1200 to 2200 m and representing different habitats, to document changes in butterfly community compositions that occurred in the past 40 years. Butterfly species are prone to fluctuate in individual numbers and our surveys were limited to three points in time (1978,2009 and 2019), but the statistical approach, the spatial and temporal scale and the focus on species occurrence (instead of species abundance) that we chose, still allow us to get crucial insights on the effects of climate and vegetation changes on butterfly community compositions.

Study Area
Butterfly communities were sampled in the SW Italian Alps (Maritime Alps, Valasco Valley) at seven sites (hereafter "plots"), scattered along an altitudinal gradient (1300-1900 m) and covering different habitat types, representative of the natural variability of the mountain system under study (Table 1).
The study area is well-known for being crucial for the conservation of the Italian insect fauna and for its rich butterfly assemblage [61]. It falls within the Site of Community Importance (SCI: "Argentera" IT1110053) included in the Maritime Alps Natural Park and it is listed among the 32 Italian priority sites for conservation known as "Prime Butterfly Areas" [62]. The position of each plot was georeferenced in 1978 [63], allowing the repetition of the surveys in 2009 and 2019 ( Figure S1). Ente di Gestione delle Aree Protette delle Alpi Marittime granted all permits necessary to perform our study within the Natural Park.

Climate Data
Historical climate data were obtained originally from the NWIOI dataset (North West Italy Optimal Interpolation-OI; [64][65][66]), processed by ARPA Piedmont (Regional Agency for the Environmental Protection), for which the ERA-40 re-analysis of the European Centre for Medium Range Weather Forecasts and raw meteorological station data were used as a background field. These data cover the period from 1958 to 2019 and include daily maximum and minimum temperature and daily cumulative precipitation. The NWIOI temperature and precipitation data were downscaled by linear interpolation to reach a resolution of 25 m, estimated to be a relevant scale for butterfly ecology. The interpolated temperature data were also corrected for elevation anomalies, using a monthly-averaged temperature lapse-rate [67]. The interpolated precipitation and elevation-corrected temperature data were finally aggregated to obtain monthly means, used in subsequent analyses. Trends and anomalies were evaluated by averaging data at annual and seasonal scales.
To assess possible altitudinal shifts of butterfly species, the climate data were used to identify two additional plots at higher elevations where vegetation and climatic conditions reflected those found at lower altitude sites in 1978. An altitudinal band representing the climatic conditions of the first year of sampling (1978) was obtained by averaging temperatures recorded in the preceding ten years (1968)(1969)(1970)(1971)(1972)(1973)(1974)(1975)(1976)(1977). The same approach was used to identify the equivalent climatic altitude band for the last decade of sampling (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018), generally at higher elevations. In this way, a new altitude band reflecting the temperature conditions found at the original altitude band of 1978 was identified, and two new additional plots were positioned in this higher-elevation band.
To analyse the time series of meteorological variables, we used two approaches. First, we adopted a Seasonal-Trend decomposition procedure based on the LOESS smoother (STL, see Cleveland et al. [68]). This method represents a filtering procedure for decomposing a time series into three components, i.e., seasonal-the periodic component present in our signal; trend-over the years; residuals-the remaining unexplained components and therefore the anomalies [68]. Second, we estimated long-term linear trends in deseasonalised temperature and precipitation time series and determined their significance using an F-test. We carried out these analyses for the full time series (1958-2019) and for two separate time periods (t1 = 1958-1979; t2 = 1980-2019), to explore possible changes in the trend intensities.
To describe the climatic characteristics of the three sampling years (1978,2009 and 2019) quantitatively, we used variables known to provide indices of referenced limitations to butterfly growth and survival rates: (i) annual daily temperature sum above 5 • C degrees from January until February, April, June and August (GDD5; surrogate for the developmental threshold for butterfly larvae); (ii) mean temperature of the coldest month (MTCO; related to overwintering survival), as suggested by Hill et al. [69] and following literature [70,71]; (iii) mean annual temperature (T mean, • C) and (iv) mean temperature for each season (DJF, winter; MAM, spring; JJA, summer).
Statistical analyses were carried out on CDO 1.9.8 (Climate Data Operator, Max Planck Institute for Meteorology 2019).

Vegetation Data
Changes in land cover were assessed through the analysis of: i) the aerial photographs taken in 1978 (source IPLA archives, Piedmont) and ii) ortho-rectified aerial images from 2010 and 2018 (source Piedmont Region administration, freely available from https:// www.geoportale.piemonte.it (accessed on 30 October 2019)). Pictures taken in 1978 were rectified on the basis of 10 ground control points retrieved from a recent (2012) orthoimage. Aerial images taken in 1978, 2010 and 2018, representative of the three periods of butterfly community sampling, were analysed for circles with a radius of 200 m around the centre of transects. Five types of land cover were identified: woodland, sparse woodland (i.e., low density of trees), scrubs, grasslands and screes. Each land cover type was quantified as a percentage of the total circle under analysis, with the exception of parts not visible in one or more images or occupied by the road. Field observations performed in 2009 and 2019, during the butterfly sampling period, were used to descriptively confirm patterns observed in the orthophotographs taken in 2010 and 2018. We compared percentages of land cover between the years and for each category by Friedman test and pairwise Wilcoxon signed-rank test. For land cover types that significantly differed between the years, we quantified the amount of change for each plot by applying Euclidean distances. First, we estimated changes for each land cover category and then for all the categories grouped, in order to obtain an aggregated measure of vegetation change per plot. To test the hypothesis that a greater change in vegetation cover determined a higher turnover in butterfly community composition, and to identify which land cover type was the most influential, we calculated Spearman rank correlation between Euclidean distances and Bray-Curtis distance matrices applied to butterfly community composition.

Butterfly Communities
Butterflies were sampled by linear transects [72], one per habitat type (plot- Figure S1), fixed in length (about 300 m) and time (45 min) in sunny days with scarce wind, between 10 am and 3 pm. All butterflies were caught by a 40 cm diameter net; only individuals hard to be identified in the field (e.g., some Pyrgus or lycaenid species) were sacrificed and identified in the laboratory. In 1978, the sampling period was limited to two sampling sessions between the end of July and the beginning of August, the ideal period to study butterflies in mountain ecosystems, because it represents the flight peak of mountain communities (Table S2). In 2009 and 2019, the sampling period was extended to cover the whole flight season from the beginning of June until the end of August, with 7 sampling sessions carried out every 7-10 days, depending on weather conditions (Table S2). The reason for lengthening our sampling time was to cover almost all the potential flight season, as to make sure that the 'no-more-found species' had really disappeared during our 40-year time frame and had not just suffered some phenological shift. The number of individuals per each species was counted for 2009 and 2019 sampling events only and used for short-time comparisons (see Sections 2.4.2 and 3.3.2). Long-time analyses (see Sections 2.4.2 and 3.3.1) are based on species occurrence; therefore, we could assess species loss or phenological shifts in a 40-year time frame but did not account for variation in species abundances.

Changes in Butterfly Communities
Our data were gathered in three sampling events far away in time (1978, 2009 and 2019). In order to avoid any bias due to natural fluctuations in the number of individuals, we based our comparison primarily on species occurrence, and we assessed community composition variations. We analysed butterfly communities sampled at Valasco valley, focusing both on the long-(LTP-1978-2019) and the short-time period (STP-2009(STP- -2019. To test LTP changes in species richness, we used a subset of the sampling sessions carried out in 2009 and 2019, which were limited to those corresponding to the time frame monitored in 1978 (two sampling sessions: 30th July and 8th August for 2009, 30th July and 10th August for 2019). Short-term variations that occurred in the last ten years were assessed by comparing the data recorded in 2009 and 2019 for the whole flight season (from the beginning of June to the end of August). We compared changes in species richness, number of plots occupied by each species and the proportion of species and of occupied plots for each ecological category among years. These analyses were carried out using ANOVA/Friedman test for LTP, and Student t-test/Wilcoxon signed-ranked test for STP according to data distribution. Shapiro-Wilk tests were used to assess data distribution. All analyses were performed using paired tests; p values were adjusted according to Benjamini-Hochberg corrections. For the STP, we assessed the variation of Shannon indexes calculated for each plot between 2009 and 2019 by applying the Student t-test for paired samples.
For the 2009 and 2019 datasets, we applied a recently developed rarity index, rr index [74] to evaluate species composition differences by detecting possible rarefaction of certain species. This index is applied to the whole study area and, following Rabinowitz's [75] scheme, it evaluates the geographical range of each sampled species by the coordinates of the observation (gri), the maximum number of the species habitat (hsi) and the maximum population size sampled (psi); the final rr index is the average of all these indices. Therefore, it functions as a local index and the classification as "rare" for a species must be intended only for the area at which it is applied (i.e., it does not mean that the species is globally or regionally rare, but it is rare in the sampled community). All of these measures range from 0 to 1, the closer the index value gets to 0, the more common the species is, while the closer it gets to 1, the rarer is the species. Following Maciel [74], rare species have an rr index higher than the average rr index calculated for the whole community. Differences in the proportion of rare species between the butterfly communities sampled in 2009 and in 2019 were tested by Wilcoxon tests.
We analysed community compositions by testing for changes in location (significant changes in community composition per plot over time) and dispersion over the years (significant changes in observed differences in community composition among plots, over time). Changes in location were tested by applying non-parametric MANOVA to Bray-Curtis distance matrices, to test if the multivariate centroids of species composition were, or were not, similar in the three groups [76,77]. Non-parametric MANOVA was performed by the function adonis of the vegan package [78]. The significance of the test was assessed by using F-tests based on sequential sums of squares obtained from permutations of the raw data (999 permutations). To account for the temporal structure and the spatial dependencies of our sampling design (7 plots at 3 points in time), we applied restricted randomisation, which did not allow for permutations across samples. Changes in dispersion were tested by the betadisper function of the vegan package, a multivariate analogous of the Levene's test for comparing group variances [76]. Non-Euclidean distances between objects and group centroids were handled by reducing the original distances to principal coordinates. To test for significance, we applied a similar randomisation approach, as previously explained.
To test for species responses over time, we used the 'indicator species' analysis (IndVal method) proposed by Dufrêne and Legendre [79]. The IndVal method combines the "specificity" of a species (its uniqueness in one of two years) and its "fidelity" (its frequency within one year). For each species, IndVal may range from 0 (no indication) to 1 (maximum indication). Statistical significance of IndVal was tested by using a Monte Carlo test, based on 999 randomisations, and performed by indicspecies package [80]. Indicator Species Analysis (IndVal method) allowed us to identify species that: To assess which species had completely disappeared over time descriptively, we used all data from the sampling sessions performed in 2009 and 2019 and compared them with data collected in 1978.

Community Temperature Index
We used the "Species Temperature Index-STI", to relate community composition to the species climatic niches [81]. The data were obtained by the updated version of the database on Lepidoptera Papilionoidea [82] build up for the project 'CkMap' by the Ministry of the Environment, Land and Sea [83]. Nowadays, this dataset includes over 388,000 individual records (the 2007 version included 60,000 records) mapped on 10 × 10 UTM grids. Temperature data were obtained by the maps of Metz et al. [84], using annual mean temperature (BIO1). The "Community Temperature Index-CTI" was calculated by averaging STI for all the species counted in each plot, for each year (see Cerrato et al. [60] for a detailed description of the procedure). CTI was used to evaluate species changes over time by applying an ANOVA/Friedman test for paired samples.
Where not otherwise specified, results have been reported in the main text as means ± standard errors. All the statistical analyses were carried out on R 3.6.1 (R Core Team 2019).  No significant trend was observed for minimum temperatures in the NWIOI data set for the study area, either for the whole time series (R 2 = 0.01, F = 1.96, df = 1, p = 0.16) or for the series divided in the two time periods (t1, R 2 = 0.01, F = 1.23, df = 1, p = 0.28; t2, R 2 = 0.06, F = 3.52, df = 1, p = 0.06). In the first years, however, minimum temperatures were usually measured at a fixed time during the night and could reflect the true minimum daily temperature only approximately. For precipitation, no trend was observed for the total period (R 2 = 0.003, F = 0.23, df = 1, p = 0.63) or for the first subset (t1 1958-1979, R 2 = 0.04, F = 0.88, df = 1, p = 0.35), in keeping with the results of Ciccarelli et al. [5]. Slight evidence of a precipitation increase was highlighted in the second period (1980-2019, R 2 = 0.24, F = 13.69, df = 1, p < 0.0001) but the R 2 value is too low to identify a real increasing trend. Comparing climatic variables between 1978, 2009 and 2019, we observed that 2009 and 2019 were hotter than 1978, but not different from each other in terms of annual mean temperatures. At a seasonal level, 2009 was a peculiar year with very low temperatures in winter, even lower than those in 1978 and much higher temperatures in spring and summer, while 2019 had a mild winter, a colder spring and a much hotter summer ( Table 2). All data are expressed as Celsius degrees ( • C). T mean = mean annual temperature; T mean (DJF) = mean winter temperature (December-January-February); T mean (MAM) = mean spring temperature (March-April-May); T mean (JJA) = mean summer temperature (June-July-August); MTCO (Jan) = Minimum Temperature of the Coldest Month, January; dd = accumulated growing degree days with a base temperature of 5 • Cfrom January until February (Feb), April (Apr), June (Jun), and August (Aug).

Vegetational Data
The comparison of the adapted aerial photographs (1978) and the orthophotographs (2010 and 2018) highlighted that significant changes in land cover occurred over time ( Figure 2). Over the period 1978-2018 we observed an impact of small wood cuttings and avalanche damage to the woodland cover, but we also recorded a significant increase in tree cover (Friedman χ 2 = 7, df = 2, p = 0.030). Between 1978 and 2010, the mean wooded surface doubled (Wilcoxon 1978-2010 z = 1, df = 2, p = 0.047), while changes between 2010 and 2018 were negligible. Fully wooded + sparsely wooded surfaces rose from 18.1 ± 4.1% to 37.7 ± 5.9%. The only observed exception was at site rh2, strongly affected by geomorphic disturbance, where the total wood cover decreased by 5.7% in the long-term period . Expanding trees were mainly represented by beech and larch at low elevations and by larch in all other sites.
An increase in the individual sizes of previously existing trees, generally determining the closure of the canopy layer, was accompanied by the presence of new tree clusters in grasslands and heathlands. In parallel to the observed increase in tree coverage, we also noticed a significant decrease (Friedman χ 2 = 11.2, df = 2, p = 0.003) in the areas occupied by grassland (1978-2010: −13.1 ± 0.1%; 1978-2018: −13.2 ± 0.1%). Scrubs also slightly decreased since 1978 (Friedman χ 2 = 4.8, df = 2, p = 0.09; Figure 2). During the last decade (2010-2018), no significant differences were assessed in the vegetation cover, the main change being linked to the density of trees, with an increase of fully wooded surfaces and a decrease of the sparsely wooded areas, whereas the sum of the two changed only slightly (+1.5%). We observed a positive and high correlation only between changes in wooded areas (Euclidean distances) and changes in butterfly communities (Bray-Curtis distances) (ρ = 0.785, p = 0.048), indicating that the expansion of wood cover enhanced dissimilarities between butterfly assemblages. No significant correlations were found between butterfly communities and other land cover types (grassland, scrubs and screes).

Changes in Butterfly Communities
The whole butterfly community recorded in the three years (considering only the two central sampling events) encompasses 77 species. In detail, 51 Table 3.

Short-Time Analysis (2009-2019)
When we compared the data collected in the whole sampling period between 2009 and 2019, we found that species richness per plot was higher in 2009 (S 2009 = 35.7 ± 3.1; S 2019 = 31.6 ± 2), but differences were not significant (t-test t = 1.619, df = 6, p = 0.156).  Table S3. Looking at rarity, in 2009 we detected 28 species (37% of total species) with rr index above the average (rr = 0.63 measured for all 76 species), while 30 species (46% of total, n = 65) were above the rr average (rr = 0.69) in 2019, identifying an increase in the proportion of rare species. The rr indices calculated for the butterfly community in 2009 and 2019 showed a slightly significant difference (Wilcoxon test W = 2927, p = 0.05). The psi index, which refers to population size, significantly decreased in 2019 (Wilcoxon test W = 3149.5, p = 0.004).

Large Scale Drivers
Ecosystems facing the coldest temperature conditions are most vulnerable to rapid climate change. These include, among others, alpine ecosystems [85]. Here, climate acts as a driver for organism distribution, altering species assemblages. Recent analyses underlined how even a moderate increase in temperature (1 • C) could lead to significant changes in arthropods community composition and species richness [49]. Our analysis of changes observed in the butterfly communities of the SW Italian Alps between 1978, 2009 and 2019 suggests the existence of scale-dependent patterns. We observed a general trend of communities changing according to large scale effects (i.e., climate change). At the same time, small-scale effects (i.e., local vegetational changes and geomorphic processes) can also potentially influence butterfly communities, acting synergistically or constituting a buffer effect.
In the study area, we observed a dramatic increase in temperatures. Using data from Open Data ARPA NWOI dataset, we identified a positive linear trend in the maximum temperatures, which could be quantified in a total increase of 2 • C during the last 40 years (+0.05 • C per year). In contrast, no clear trend was detected in precipitation patterns, and minimum temperatures practically did not vary. Such observed patterns correspond rather well to the trends in climatic variables summarised by Ciccarelli [5] and by Acquaotta [6], although these works observed shallower increases in maximum temperatures (respectively +0.023 • C and +0.018 • C per year) compared to our findings. A reason for this may be that these authors analysed the whole western Alpine region, encompassing wider altitudinal and latitudinal ranges. Instead, we focused on smaller and finer scale, at a level relevant for butterfly communities. Still, since even a small increase in maximum temperatures has potential implications on the composition of butterfly assemblages, by changing the distributional limits of some species, or by promoting the expansion or the local extinction of others [86], it is not surprising that in our study area this temperature increase impacted the biodiversity of the whole butterfly community, mainly in terms of species composition according to their functional traits.

Small Scale Drivers
A generalised spread of forested areas (open woodland and dense woodland) and a reduction in the size of grasslands during the past decades has been observed for all our study areas. Heathland expansion on the previous grasslands was paralleled by the closure of the canopy layer of the larch forest, accompanied by changes in dominant species. In particular, we observed an increase in some small and tall shrubs species, which were by far less represented in the few qualitative observations from 1978 (e.g., Juniperus nana, Laburnum anagyroides) and a decrease in abundance of Rhododendron ferrugineum. The areas less subjected to the forest cover increase are located near the last remaining centre of pastoral activity (rh2, wm, pf). In all other transect areas, we observed a marked growth of woodland occupancy, with the exception of ss, due to the absence of soil in the block scree.
In the whole study area, and especially in two of the transects (rh1, rh2), we detected cyclic habitat modifications due to periodic avalanches. In the last decade, we observed a high abundance of Rubus idaeus in these sites, mainly represented by one-year stems only. The spread and survival of raspberry in the study areas was certainly facilitated by geomorphic disturbance, as known for opportunist species of the subalpine and alpine communities [23,87]. In the European Alps, the effect of climate change, a large-scale driver, is confounded by local-scale human activities. Cattle grazing in the alpine pastures has decreased throughout the last century, allowing a fast recolonisation by trees and shrubs, but this can also result from climatologic variation [21,25,54,88].

Changes in Butterfly Assemblages
The modifications observed in the composition of species assemblages apparently go in the direction of the predicted impacts of temperature increase. In particular, we wish to stress the loss or strong decline in a number of species whose characteristics make them potentially sensitive to climate change, e.g., the geographically localised, the poor dispersers or the ecologically highly specialised butterflies [89,90]. At the same time, we also reported the expansion of some species, such as woodland, thermophilous and generalist, as already recorded at the European scale and generally associated with a global increase in temperature [14,46,54,91,92]. More in detail, since 1978 we observed: (i) a continuous increase in maximum temperatures ( Figure 1); (ii) an increase of forest coverage paralleled by a decrease of grasslands ( Figure 2); (iii) changes in community compositions ( Figure 6) determined by non-random turnover.
Communities sampled in 1978 significantly differed from those sampled in 2009 and in 2019. We observed the loss of specialised species (i.e., Coenonympha glycerion, Colias phicomone, Pieris callidice). At the same time a higher number of the truly alpine species (i.e., those only flying above or at the tree line, such as Colias phicomone or Pieris bryoniae) were only observed in 1978. Conversely, the last decade (2009-2019) shows an increase in relatively thermophilous elements (i.e., Satyrus ferula, Parnassius apollo) and in shrubland/woodland species (i.e., Argynnis paphia, Brenthis daphne). Among the species protected at European level by the Habitats Directive (92/43/EEC), P. apollo and P. mnemosyne showed changes in occupancy and flight period, respectively. Both these two charismatic species are suffering decline in other parts of Italy [93]. Significant results of non-parametric MANOVA (Figure 6a,b) show that butterfly communities have changed across years and that species compositions recorded in 2019 and 2009 were not analogous to those observed in 1978, but represent a reassembly of the original species along with the ingression of several new elements (only 23 species to 77 were found in all years). Indeed, many species showed no stable trend through time, due to changes in occupancy, elevation or phenology (Figures 3 and 4). Nevertheless, despite the continuous increase of temperatures, the comparison between communities sampled in 2009 and 2019 shows minor differences (Figure 6c) in the position of their centroids. This seeming local resilience is mainly due to the inflow of 11 'new' and relatively generalist species in the short period (2009-2019). More in general, we noticed a tendency towards biotic homogenisation referred to the increase in biological similarity among communities (sensu [94]). Group variances between years were not significantly different (see the multivariate analogous of the Levene's test in 3.3.1 result section). However, a reduction in the distinctiveness of sampled communities has been observed, in particular from 1978 to 2009 ( Figure 6).
Even if studies on the community-level responses of butterflies through time are less frequent than analyses of single species responses, results similar to ours are paving common routes and patterns in butterfly ecology. For example, a similar change in community composition over time, along with an increase in community similarity, has been observed in the data analysis from the UK Butterfly Monitoring Scheme over a 20-year period [95]. In this latter case, the main causes were related to an increase in southerly distributed species, particularly the more generalist ones, which spread northwards, accompanied by the decline of some localised species. In mountain ecosystems, Wilson et al. [92] observed that widespread species became increasingly dominant during the last 40 years, while highaltitude species strongly declined. Similar patterns were observed along the altitudinal gradient in the space-for-time studies even on a short time gap [60]. Climate change acts by filtering vulnerable species as observed for habitat simplification [96], but this phenomenon was already observed as a consequence of habitat deterioration, such as in the case of agricultural intensification [97]. This rearrangement at the community level is accompanied by some small changes in the relative representation of some ecological groups. Looking at the ecological categories of the sampled species, we can observe that differences are mainly related to an increase in the proportion of low-altitude species and changes related to habitat specialisation. Species linked to the wooded habitats strongly increased in frequency, apparently at the expenses of species linked to grasslands occurring above the tree line (subalpine and alpine belts). We can also observe a generalised ingression of widespread and common species (e.g., Ochlodes sylvanus, Aglais urticae, Papilio machaon, Colias crocea, Pieris brassicae; see Figure 3 and Table S1), mainly characterised by high vagility and broad ecological tolerance, as well as a general spread of thermophilous species, as confirmed by an increase in the Community Temperature Index ( Figure 5).
More interestingly, the analysis of indicator species by the IndVal method (Table 3) and the comparison of the list of species observed in 1978 with the full dataset of 2009 ( Figure 3, Table S1) identifies some 'winner' species (the expanding ubiquitous species in our dataset), as opposed to others that may be viewed as 'losers' (replaced species) [98]. In particular two species were completely absent from the 1978 dataset and were present in most of the plots sampled in 2009: Argynnis paphia (in 5 of 7 plots) and Brenthis daphne (in all plots). These species are characterised by a long flight period, spanning from the beginning of July until the middle of August, and would certainly have been recorded, if at all present, in 1978. Both species have undergone, during the past decades, a northwards expansion accompanied by a stable southern boundary [91,99,100], while an increase in the elevational optimum has also been recorded in the case of Argynnis paphia [8]. Argynnis paphia is a woodland butterfly, which uses as larval host plant some Viola species, but lays its eggs in the tree bark, where the newly hatched larvae hibernate, and selects sites where wellshaded larval host plants are abundant. Brenthis daphne is a sub-nemoral, scrubland species and uses as larval host plants many species of Rubus. Brambles are early successional, highly competitive, fast-growing species that are expanding everywhere in Europe, being generally favoured by an abundant nitrogen supply; in mountain ecosystems, they can quickly occupy places set free from avalanches [101,102].
To come now to the 'losers', we observed three cases of apparent species substitution (indicating with this term the complete disappearance of one species, replaced by a congeneric one since 1978) at our sampling sites, indicating a general loss of specialised and narrow-range species and a general expansion of ecologically tolerant elements. In the last decade, we neither recorded Coenonympha darwiniana, an alpine endemic, nor Coenonympha glycerion, a hygrophilous species, but we found Coenonympha arcania, that had not been recorded in 1978. The latter is a widespread species, characteristic of the lowland ecotonal habitats and a much stronger generalist than the previous two. A similar case is that of Colias phicomone, a xerophilous species, linked to open herbaceous habitats occurring above the tree line (subalpine and alpine belts). It was repeatedly found in 1978 but was substituted in 2009 by the ingression of Colias crocea, an ecotonal, generalist and rather thermophilous species, characterised by very high vagility. Finally, Pieris callidice, a microthermic, xerophilous species of the open herbaceous environments of the alpine belt, has disappeared completely from the study sites, where we found Pieris daplidice, a generalist, thermophilous species, typical of the low elevations. Colias phicomone and Pieris callidice did not show any altitudinal shift since they have not been found in higher elevations plots sampled in 2019.
As a word of caveat, we note that one source of bias in our analysis is represented by the reduced sampling period in 1978. The strong observed differences in community composition might be reduced if sampling data from 1978 were available for all the favourable season (June-August), allowing to disentangle real differences in community composition from simpler changes in the flight period of individual species. To overcome this, we selected a much longer time frame for samplings carried out in 2009 and 2019, so that species having anticipated their flight period would have been detected anyway, but of course we could not use the same procedure for 1978. In the latter years, some species might have flown before, or perhaps later, than our sampling period. This is, however, unlikely. On the one hand, our climatic analysis showed that 1978 was colder than 2009, and species would not tend to fly earlier in such conditions. On the other hand, early to mid-August is the end of the flying period of butterflies throughout the Italian Alps, with very few exceptions. In any case, such a reduced sampling period (end of July-beginning of August) represents the peak of activity for butterfly communities in the mountain ecosystems of the Alps. Moreover, the prolonged sampling period in 2009 and 2019, which covered the whole favourable season, has allowed to reveal that some (at least 8) species have completely disappeared. Secondly, we observed that most species present in the reduced dataset of 2009 and 2019 (end of July-beginning of August, used for the statistical comparisons) and missing from the 1978 dataset are species characterised by long flight periods, centred in the middle of summer. Consequently, they should have been found, if present, also in 1978. Such observations confirmed the robustness of our results and the reliability of changes in community compositions through time.

Conclusions
Our work represents the first attempt to describe changes in butterfly community composition in the European Alps over a long-term temporal scale. Responses to climate and habitat changes vary widely among species occurring in the same communities, and consequently it is important to understand the heterogeneity of species responses and its implication on changes in community compositions. However, it is difficult to assess the relative role of climate and vegetation modifications in determining the observed pattern in butterfly communities. The observed sharp increase in vagile, tolerant and thermophilous species perfectly fits the consequence of an increase in temperature ( [41,57,92,103,104] and references therein). The reduction of open areas, in fact, is not a problem for the highly mobile species that can easily reach another suitable site, while the less mobile species will remain locked to their continuously shrinking micro-environment [105,106]. In the same way, because of their strong adaptability, the broadly tolerant species will take advantage of a changing environment, to the detriment of the stenotopic and more specialised species, thereby contributing to reduce the community distinctiveness [106,107].
Species respond individually to the changing environmental conditions, depending mainly on their functional traits and habitat requirements [108,109]. This can originate new species assemblages, which can be appreciated only by the examination of entire communities throughout time (e.g., [92,110,111]). We observed the highest changes during the first time frame , while in the last decade a slow and background erosion of species emerged. In this framework, we can envisage that long-term environmental change can enhance the observed pattern of alterations in butterfly guild compositions consistently with predicted impact by local warming.
Several studies investigating community changes over time relied on comparisons between current data and historical data sets (atlases, collection specimens), the latter of which were generally collected in a non-standardised way and/or referred to a much coarser spatial grain [112]. Transects set in well-specified areas represent a more appropriate tool, e.g., [60]. The implementation of projects aiming to record and share data on butterfly observations, such as the Italian Butterfly Monitoring Scheme (https://butterfly-monitoring.net/it (accessed on 23 November 2021), will guarantee the availability of robust and standardised data for research addressing long-or short-term butterfly assemblage variations.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/insects13010043/s1, Table S1: Butterfly community sampled in 1978, 2009 and 2019;  Table S2: Information about sampling events; Table S3: Indicator species (IndVal species) for both communities; Figure S1: Location and altitudinal gradient of sites surveyed by transects (pollard walk) in the Valasco Valley (SW Italian Alps); Figure S2: Frequency distribution of the species-by-species number of "gained" or "lost" plots.